Rozważmy problem Stokesa, który wymaga dodatkowej stabilizacji. Zapiszmy go w postaci pojedynczego funkcjonału
\( B(({\bf u},p),({\bf v},q)) = L({\bf v }) \)
gdzie
\( B(({\bf u},p),({\bf v},q)) = \left( \begin{matrix} a({\bf u}, {\bf v}) + b(p, {\bf v}) \\ b(q, {\bf u} ) \\ \end{matrix} \right) \)
\( L({\bf v }) = \left( \begin{matrix} f({\bf v} \\ 0 \\ \end{matrix} \right) \)
oraz
\( a({\bf u}, {\bf v }) =\int_{\Omega} \frac{\partial u_1}{\partial x}\frac{\partial v_1}{\partial x}dxdy + \int_{\Omega} \frac{\partial u_1}{\partial y}\frac{\partial v_1}{\partial y}dxdy + \int_{\Omega} \frac{\partial u_2}{\partial x}\frac{\partial v_2}{\partial x}dxdy + \int_{\Omega} \frac{\partial u_2}{\partial y}\frac{\partial v_2}{\partial y }dxdy \)
\( b(p, {\bf v }) = \int_{\Omega}p \frac{\partial v_1}{\partial x} dxdy+ \int_{\Omega} p \frac{\partial v_2}{\partial y } dxdy \)
Oznaczamy przez \( {\bf U} \times P \) przestrzeń funkcji w którym poszukujemy rozwiązanie problemu Stokesa \( ({\bf u},p) \in {\bf U } \times P \).
Następnie oznaczamy przez \( {\bf V } \times Q \) przestrzeń funkcji których używamy do testowania w naszym sformułowaniu wariacyjnym (zwanym też słabym).
W rozdziale tym opiszemy za pomocą trochę abstrakcyjnego języka matematycznego jak wyprowadza się metodę minimalizacji reziduum [1].
Nasz oryginalny problem zapisać można tak: szukamy funkcji \( ({\bf u},p) \in {\bf U } \times P \) takich że \( B( ({\bf u},p),({\bf v},q) ) = L({\bf v }) \) dla wszystkich funkcji testujących \( ({\bf v },q) \in V\times Q \).
Matematycznie można oznaczyć lewą stronę \( B(u,p) \) a prawą stronę \( L \) wprowadzając pewne operatory matematyczne (matematycy powiedzą że \( B : ({\bf U}\times P) \rightarrow (V\times Q)' \) oraz \( \langle B({\bf u},p), ({\bf v},q) \rangle =B(({\bf u},p),({\bf v },q)) \).
Wówczas nasz problem sprowadza się do znalezienia \( ({\bf u},p) \) takiego że
\( B({\bf u },p)-l=0 \). Wyrażenie to nazywa się reziduum \( r=B({\bf u },p)-l \).
Jeśli mamy jedynie rozwiązanie przybliżone, na przykład oznaczone \( ({\bf u_h },p_h) \)
wówczas nasze reziduum będzie wyrażało błąd, im większe będzie reziduum tym błąd rozwiązania przybliżonego będzie większy \( r_h=B({\bf u_h },p_h)-l \).W jaki sposób zmierzyć błąd rozwiązania, czyli wielkość reziduum?
Żeby zmierzyć jak duże jest reziduum musimy wprowadzić jakąś normę (czyli miarę wielkości) \( r_h = \| B({\bf u_h },p_h)-l \| \). Mogę teraz powiedzieć że szukamy takiego rozwiązania naszego problemu \( ({\bf u_h },p_h) \) żeby reziduum (czyli błąd) był jak najmniejszy. Matematycznie zapisuje się ten warunek w postaci minimalizacji normy
\( ({\bf u_h },p_h) = \displaystyle{argmin}_{({\bf w_h },s_h) \in U\times P } {1 \over 2} \| B({\bf w_h},s_h) - l \|_{({\bf V},Q)' }^2 \)
Przed normą dopisuje się \( 1\over 2 \) oraz podnosi się normę do kwadratu.
Problem praktyczny jest taki że nie wiemy jak policzyć normę z różnicy \( \|B({\bf w_h },s_h)-l\| \). Normę tą mierzy bowiem odległości w abstrakcyjnej przestrzeni matematycznej \( (V\times Q)' \) (matematycy powiedzą że przestrzeń ta jest przestrzenią dualną do przestrzeni funkcji testujących).
Żeby rozwiązać ten problem, wykonujemy trick matematyczny polegający na rzutowaniu normy z przestrzeni \( (V\times Q)' \) do przestrzeni testującej \( V\times Q \). Wprowadza się operator przerzucający przestrzeń \( (V\times Q) \) w przestrzeń \( (V\times Q) ' \). Operator ten nazywamy operatorem Riesza \( R_{(V,Q)} \colon (V\times Q) \ni (v,q) \rightarrow ( ( v,q),.) \in (V\times Q)' \). Operator odwrotny natomiast przerzuca przestrzeń \( (V\times Q)' \) z powrotem w przestrzeń \( V\times Q \), czyli \( R_{(V,Q)}^{-1} \colon V' \ni (v,.) \rightarrow v \in V \). Tak więc jeśli użyjemy odwrotnego operatora Riesza, przerzucimy nasz problem do przestrzeni funkcji testujących \( V\times Q \). Mamy więc
\( ({\bf u_h},p_h) = \displaystyle{argmin}_{({\bf w_h },s_h) \in (U_h\times P_h)} {1 \over 2} \| {R_{(V,Q)}^{-1}} (B({\bf w_h},s_h) - l) \|_{V\times P }^2 \). Szukamy więc punktu w którym funkcja osiąga minimum. Taki punkt osiągany jest wtedy kiedy pochodne tej funkcji są równe zero. Ponieważ mamy tutaj funkcjonał, matematycy powiedzą że szukamy takiego \( ({\bf u_h },p_h) \) w którym pochodne Gateaux są równe zero we wszystkich kierunkach
\( \langle R_{(V,Q) }^{-1} (B({\bf u_h },p_h) - l), R_{(V,Q)}^{-1 }(B,({\bf w_h },s_h) ) \rangle_V = 0 \quad \forall ({\bf w_h },s_h) \in U_h \times P_h \) gdzie \( ({\bf w_h},s_h) \) oznacza wszystkie możliwe "kierunki" ze skończenie wymiarowej przestrzeni funkcji \( U_h \times P_h \) (w której przy okazji szukamy rozwiązania naszego problemu Stokesa).
Formalnie rzecz ujmując, w tym miejscu musimy skonstruować skończenie wymiarową przestrzeń funkcji \( U_h\times P_h \) w której będziemy szukać naszego rozwiązania. Robimy to oczywiście w naszym podręczniku stosując funkcje bazowe B-spline, rozpinając je zgodnie z wektorami węzłów wzdłuż osi x oraz wzdłuż osi y.
Wzdłuż osi \( x \) wprowadzamy wektor węzłów [0 0 0 1 2 3 4 4 4], podobnie wzdłuż osi \( y \) wprowadzamy wektor węzłów [0 0 0 1 2 3 4 4 4].
Uzyskaliśmy dwie bazy jednowymiarowych funkcji B-spline \( \{ B_{i,2}(x) \}_{i=1,...,6 } \) oraz \( \{B_{j,2}(y)\}_{j=1,...,6 } \). Następnie utworzymy z nich iloczyny tensorowe \( B_{i,j;2}(x,y)=B_{i,2}(x)B_{j,2}(y),i,j=1,...,6 \).
Zauważmy że nasz obszar \( \Omega \) rozpięty jest na kwadracie \( [0,1]^2 \), podobnie jak nasze 6*6=36 funkcji bazowych co wynika z definicji wektora węzłów [0 0 0 1 2 2 3 4 4 4].
Nasza przestrzeń aproksymacyjna jest więc rozpięta przez iloczyny tensorowe naszych jednowymiarowych funkcji B-spline. Teraz mamy natomiast dwie współrzędne prędkości oraz jedną wartość skalarną ciśnienia, dostajemy więc
\( U_h \times Q_h = gen \{ B_{i,j;2}(x,y)=B_{i,2}(x)B_{j,2}(y),i,j=\\=1,...,6, B_{i,j;2}(x,y)=B_{i,2}(x)B_{j,2}(y),i,j=1,...,6, B_{i,j;2}(x,y)=B_{i,2}(x)B_{j,2}(y),i,j=1,...,6 \} \).
Oznaczamy
\( r=R_{(V,Q)}^{-1}(B({\bf u_h },p_h)-l) \)
i wówczas nasz problem minimalizacji reziduum zapisać można w postaci
\( \langle r , R_{(V,Q)}^{-1} (B\, ({\bf w_h },s_h) ) \rangle = 0 \quad \forall \, (w_h,s_h ) \in U_h\times P_h \)
co jest równoznaczne
\( \langle B({\bf w_h },s_h), r \rangle = 0 \quad \quad \forall (w_h,s_h ) \in U_h \times P_h. \)
Dorzucamy teraz drugie równanie, korzystając z definicji reziduum \( r=B({\bf u_h },p_h)-l \), przemnażam definicje reziduum przez funkcje testujące z przestrzeni funkcji testujących i dostaje równanie:
\( (r,v)_V=\langle B({\bf u_h },p_h)-l,(v,q) \rangle \quad \forall( v,q )\in V\times Q. \)
Mam więc układ równań: znajdź reziduum w nieskończenie wymiarowej przestrzeni funkcji testujących \( r\in V\times Q \) oraz rozwiązania w skończenie wymiarowej przestrzeni funkcji \( {\bf u_h },p_h) \in U_h\times P_h \)
\( (r,v)_V - \langle B({\bf u_h },p_h)-l ,(v,q) \rangle = 0 \quad \forall (v,q) \in V\times Q \\ \langle B({\bf w_h },s_h),r\rangle = 0 \quad \forall (w_h,s_h) \in U_h\times P_h \)
Gdybyśmy byli wstanie rozwiązać to równanie, to dostali byśmy najlepszą możliwą aproksymacje naszego problemu w przestrzei funkcji aproksymujących \( U_h\times P_h \). Niestety nie jest to możliwe, ponieważ nieskończona przestrzeń funkcji testujących generuje nam nieskończenie wiele równań które należało by rozwiązać. Musimy więc wybrać drugą inną bazę do przybliżenia przestrzeni testowej. Tutaj pojawia się pierwsza istotna różnica metody minimalizacji reziduum z tradycyjną metodą elementów skończonych. Otóz w metodzie minimalizacji reziduum mamy dwie różne przestrzenie, jedną do aproksymacji rozwiązania, drugą do testowania. Oczywiście przestrzeń testującą również rozpinamy za pomocą funkcji B-spline. Na przykład możemy użyć takiego samego patcha elementów jak dla przestrzeni aproksymującej (ale nie jest to wymagane!), ale podnieść stopień funkcji B-spline.
Precyzujemy jak wyglądać będzie nasz patch elementów oraz nasze funkcje bazowe B-spline rozpięte na elementach, podając wektor węzłów wzdluż osi \( x \) oraz wektor węzłów wzdłuż osi \( y \). Odsyłamy tutaj do rozdziału trzeciego w którym opisany jest sposób definiowania funkcji bazowych za pomocą wektorów węzłów i wielomianów B-spline.
Wzdłuż osi \( x \) wprowadzamy wektor węzłów [0 0 0 1 2 3 4 4 4], podobnie wzdłuż osi \( y \) wprowadzamy wektor węzłów [0 0 0 1 2 3 4 4 4]. Używamy więc do aproksymacji przestrzeni kwadratowych B-spline'ów o ciągłości C1.
Do testowania zwiększymy ich stopień, oraz pozostawimy ciągłość C1. Wzdłuż osi \( x \) wprowadzamy wektor węzłów [0 0 0 0 1 1 2 2 3 3 4 4 4 4], podobnie wzdłuż osi \( y \) wprowadzamy wektor węzłów [0 0 0 0 1 1 2 2 3 3 4 4 4 4].
Generalnie przestrzenie testujące muszą być większe niż przestrzeń aproksymacyjna. Możemy uzyskać ten efekt zwiększając liczbę elementów lub podnosząc stopień funkcji B-spline, lub redukując ciągłość bazy funkcji B-spline, lub mieszając te trzy podejścia razem. Przestrzeń aproksymacyjna powinna być podprzestrzenią przestrzeni testującej.
Uzyskaliśmy dwie bazy jednowymiarowych funkcji B-spline stosowanych do testowania \( \{ \tilde{B} _{k,3}(x) \}_{k=1,...,10 } \) oraz \( \{\tilde{B}_{l,3}(y) \}_{l=1,...,7 } \). Następnie utworzymy z nich iloczyny tensorowe \( \tilde{B}_{k,l;3}(x,y)={\tilde{B}}_{k,3}(x){\tilde{B}}_{l,3}(y),k,l=1,...,10 \).
Nasza przestrzeń testujące, uwzględniająca dwie współrzędne funkcji testującej dla pola prędkości oraz jedną współrzedną funkcji testującej dla ciśnienia, to
\( V_h\times Q_h = gen \{ \tilde{B}_{k,3}(x)\tilde{B}_{l,3}(y),i,j=1,...,10; \tilde{B}_{k,3}(x)\tilde{B}_{l,3}(y),i,j=1,...,10; \tilde{B}_{k,3}(x)B_{l,3}(y),i,j=1,...,10 \} \).
Mając dyskretną przestrzeń testującą, uzyskujemy wreszcie nasz dyskretny (czyli skończenie wymiarowy) układ równań który będziemy chcieli rozwiązać,
Szukamy \( ( r_m, ({\bf u_h },p_h) )_{ V_m \times U_h \times P_h } \)
\( (r_m,({\bf v_m }, q_h))_{V_m \times Q_h } - \langle B ({\bf u_h },p_h)-l,({\bf v_m }, q_m) \rangle = 0 \quad \forall ({\bf v_m },q_m) \in V_m \times Q_m \\ \langle B({\bf w_h }, q_h),r_m\rangle = 0 \quad \forall ({\bf w_h },q_h) \in U_h\times P_h \)
Pojawia się tutaj tak zwany iloczyn skalarny (wewnętrzny) w dyskretnej przestrzeni testującej \( (*,*)_{V_m \times Q_M } \).
Żeby odgadnąć jaki iloczyn skalarny użyć w naszym problemie, musimy popatrzeć na nasze oryginalne równanie zapisane w formie słabej
\( B( ({\bf u}p),({\bf v}, q)) = \left( \begin{matrix} a({\bf u}, {\bf v}) + b(p, {\bf v}) \\ b(q, {\bf u } ) \\ \end{matrix} \right) \)
i zobaczyć w jakiej formie są funkcje testujące \( ({\bf v },q) \) w jej składowych
\( a({\bf u}, {\bf v }) =\int_{\Omega} \frac{\partial u_1}{\partial x}\frac{\partial v_1}{\partial x}dxdy + \int_{\Omega} \frac{\partial u_1}{\partial y}\frac{\partial v_1}{\partial y}dxdy + \int_{\Omega} \frac{\partial u_2}{\partial x}\frac{\partial v_2}{\partial x}dxdy + \int_{\Omega} \frac{\partial u_2}{\partial y}\frac{\partial v_2}{\partial y }dxdy \)
\( b(p, {\bf v }) = \int_{\Omega}p \frac{\partial v_1}{\partial x} dxdy+ \int_{\Omega} p \frac{\partial v_2}{\partial y } dxdy \)
Nasz produkt wewnętrzny powinien więc zawierać pochodne funkcji \( v_1,v_2 \) oraz wartości funkcji \( q \).
\( ( ({\bf u },p), ({\bf v },q))_{V_m\times Q_m} =\\= \int_{\Omega} p(x,y) q(x,y) dxdy + \int_{\Omega} \frac{\partial u_1(x,y)}{\partial x } \frac{\partial v_1(x,y) }{\partial x }dxdy + \int_{\Omega} \frac{\partial u_1(x,y)}{\partial y } \frac{\partial v_1(x,y)}{\partial y } + \int_{\Omega} \frac{\partial u_2(x,y)}{\partial x } \frac{\partial v_2(x,y)}{\partial x }dxdy + \int_{\Omega} \frac{\partial u_2(x,y)}{\partial y } \frac{\partial v_2(x,y)}{\partial y }dxdy \)
Podsumowując, w metodzie minimalizacji reziduum, musimy zdefiniować osobną przestrzeń aproksymacyjną, oraz osobną (większą) przestrzeń testującą, oraz wybrać iloczyn wewnętrzny przestrzeni testującej. Dostajemy wówczas układ równań w którym niewiadome to nasze rozwiązanie \( ({\bf u },p) \) oraz dodatkowo reziduum \( r \).
To jak dobrze działać będzie metoda minimalizacja reziduum zależy w dużej mierze od naszego wyboru przestrzeni testującej oraz iloczynu wewnętrznego. Jeśli iloczyn wewnętrzny będzie dostatecznie dokładny, oraz przestrzeń testująca będzie dostatecznie duża, wszystko będzie idealnie działać, i dostaniemy najlepsze możliwe do uzyskania rozwiązanie w przestrzeni aproksymacyjnej \( U_h \) (ale nie lepsze niż pozwala na to przestrzeń aproksymacyjna). Matematyczne uzasadnienie jest takie, że problem minimalizacji reziduum z nieskończenie wymiarową przestrzenią testową realizuje warunek stabilności inf-sup ze stałą równą 1. Jeśli przestrzeń testującą zawęzimy do przestrzeni skończenie wymiarowej, to warunek inf-sup może nie być już idealnie realizowany w tej przestrzeni. Zwiększanie przestrzeni testującej będzie nam poprawiać realizowalność warunku stabilności inf-sup.
Nasz układ równań do rozwiązania metodzie minimalizacji reziduum wyglada następująco
\( \left( \begin{matrix} G & B \\ B^T & 0 \\ \end{matrix} \right) \left( \begin{matrix} r \\ u \end{matrix} \right)= \left( \begin{matrix} f \\ 0 \end{matrix} \right) \)
Zobaczmy jak wygląda nasz układ równań na przykładzie konkretnych przestrzeni aproksymacyjnej i testującej.